%--------------------------------------------------------------------------
%
% tries to find an Ma such that the eigenvalue lambda = x + 1i * y, where
% the real part, x, is prescribed.


function Ma = find_complex_eigs

x = -120;

h = 0.7;
k = 1e-1;
p = params;

opts = optimset('TolFun', 1e-10);
X = fsolve(@(X) fun(X, x, h, k, p), [2e5; 5], opts);

Ma = X(1);

fprintf('Ma = %.4e, Im{lambda} = %.4e\n', X(1), X(2));



function F = fun(X, x, h, k, p)

Ma = X(1);
y = X(2);

lambda = x + 1i * y;
mu = sqrt(-k^2 - lambda);

par = p;
par.Ma = Ma;

A = w_coeffs(1, h, k, par);
w = @(z)  (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);

v_0 = (1 - p.beta) * (1 - 1 / h);
c = p.beta + v_0 - p.delta * (1-p.beta)*(p.beta + v_0) / 2;
gamma = p.delta * (1 - 2 * c);

tmp = mu * sin(mu * h) - gamma * cos(mu * h) + ...
    p.delta * (1-p.beta) * (p.beta + v_0) / h^2 * quad(@(z) w(z) .* z .* cos(mu*z), 0, h);


F(1) = real(tmp);
F(2) = imag(tmp);